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Abstract. In this paper, we discuss theoretically the behavior of the four point 
nonlinear susceptibility and its associated correlation length for supercooled liquids 
close to the Mode Coupling instability temperature T c . We work in the theoretical 
framework of the glass transition as described by mean field theory of disordered 
systems, and the hypernetted chain approximation. Our results give an interpretation 
framework for recent numerical findings on heterogeneities in supercooled liquid 
dynamics. 
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1. Introduction 

Recently, a lot of attention has been devoted to understanding the nature of dynamical 
heterogeneities in supercooled liquids Q, H [!, U, !, H ■ Many numerical experiments 
have found long lived dynamical structures which are characterized by a typical length 
and a typical relaxation time which depend on the values of the external parameters 
(temperature and density). A way to quantify this dynamical heterogeneities is in terms 
of the 4-point density function, and its associated non-linear susceptibility, which show 
power law behavior as one approaches the Mode Coupling temperature T c from above. 
In this paper we review the details of the theoretical calculations of this function put 
forward in || |9], [TOj] and discuss some new results. 



At the glass transition one observes freezing of density fluctuations. The function 

92{x) = (00 + y) - p){p{y) -p)> (1) 

is often regarded as the Edwards-Anderson order parameter signaling the onset of 
glassiness. It is therefore quite natural to try to interpret the dynamical heterogehinities 
and the correlation length in terms of fluctuations of the order parameter, and study 
the 4-point function 

g A (x) = ([(p(x +y)- p)(p(y) - p)] 2 > - ((p(x + y) - p)(p(y) - p)> 2 , (2) 
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and its related non-linear susceptibility Xi — I dxg±(x). To our knowledge the first 
proposal to study the 4-point function to identify a growing correlation length in 
structural glasses was in [ 1 1 1 in the context of a numerical study of a Lennard- Jones 
liquid. There no sign of growing correlation was found, probably because of the 
insufficient thermalization. However, more accurate measurements j|, |5], |9|, |TIJ show 
that there is a dynamical correlation length which grows as T c is approached. 

Here we would like to investigate theoretically the behavior this function in the 
context of the picture of the glass transition that comes out from the study of disordered 
mean-field models |13|], and from some approximation scheme of molecular liquids [14[ . 

In mean field disordered systems one finds that decreasing the temperature from the 
liquid phase, two different transitions appear: a dynamical transition at a temperature 
T c , and a static (Kauzmann-like) transition at a lower temperature Tk- At the 
dynamical transition T c , identified with the Mode Coupling Theory ]15| transition 
temperature, equilibrium density fluctuations freeze and ergodicity breaks down. Below 
that temperature, the Boltzmann distribution is decomposed in an exponentially large 
number of ergodic components e NTi ( T ) . S(T), the logarithm of the number of these 
components is the configurational entropy, which decreases for decreasing temperatures, 
and the "static" transition signals the point where S(Tr-) = 0. Dynamically, a non zero 
Edwards-Anderson order parameter signals freezing. 

As it has been many times remarked, this theory misses the existence of 
local activated processes which restore ergodicity below T c . These can be included 
phenomenologically to complete the picture. We will suppose that the ergodic 
components which the ideal theory predicts below T c become in real systems metastable 
states (or quasistates), capable to confine the system for some large, but finite times on 
given portions of the configuration space. The inclusion of activated processes, although 
done by hand, has far reaching consequences. 

The foundation of the notion of quasistates is based on the time scale separation 
(as it can be seen in the shape of the structure function), which allows to consider 
"fast" degrees of freedom quasi-equilibrated, before the "slow" degrees of freedom can 
move. So, this notion applies below as well as above T c , where the two step relaxation 
is predicted even by the ideal theory. This point has been recently stressed in (T^j in 
a different context. Both above and below T c we can talk of quasistates in which the 
system equilibrates almost completely before relaxing further. The typical life time of 
the quasistate will be of the order of the alpha relaxation time r a . 

Our basic observation is that within the described theoretical framework, the 
quasistates correspond to highly correlated regions of the configuration space, typical 
configurations belonging to the same quasistate would appear to be highly correlated. 
On the other hand, configurations belonging to distant quasistates as the ones which 
correspond to large time separation t » r a , show typically low correlations. 

We argue then, that the dynamical correlation length and susceptibility observed 
in the simulations referred to above, can be estimated by the corresponding quantities 
within a quasistate. On the other hand, the long time limit of the same quantities, i.e. 
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the value reached for times much larger then the lifetime of the quasistates, correspond 
to maximally distant quasistates. This predicts maximal fluctuations and heterogeneity 
on a time scale of the order of r a . 



2. How to compute quasistate averages: "recipes for metastable states" 

In this section, we address the question on how to compute correlation functions within 
singles quasistates, reviewing some "recipes" that were put forward in [|T7| , 18| . Let us 



consider the case of mean-field spin glass models below T c , where there is true ergodicity 
breaking and the quasistates are true ergodic components. Suppose to be above T K 
so that the configurational entropy S(T) > 0. Given any local observable A(x), its 
Boltzmann average can be decomposed as 

(A(x)) 

Boltzmann 

^2w a (A(x)) a , (3) 

a 

where the index a runs over all the e Nj: states, the weights of the different states 
w a would all be of the same order w a ~ exp(— iVE(T)). In the following we will be 
interested to compute space averages (correlation functions) among local observables, 
/ dx{A(x))(B(x + y)). If by (■) we mean Boltzmann average, we can expand each of the 
two averages according to (|3|) and find that 

J dx(A(x))(B(x + y)} = J dxY,w a w p (A(x)) a (B(x + y)) p , (4) 

a,/3 

which, due to the fact that the number of ergodic components is exponentially large, is 
dominated by the terms in the double sum with a ^ (3. 

Our major interest will be to compute instead averages of the kind 
/ dx J2 a w a{A(x)) a (B(x + y)) a i.e. correlation function in the same ergodic component. 



To this purpose we can use a conditional Boltzmann prescription ||I7|, |Tq| , where one fixes 
a reference configuration Y = {y%, yisr}, and only the configurations X = {x%, xn} 
similar enough to the reference configuration are given a non vanishing weight. 

Let us consider as a notion of similarity among two configurations X and Y the 
function, that we call overlap, 



Q(X,Y) = j dxdy(p x (x) - p)(p Y (y)- p)w(\x-y\) (5) 

where: 

• pz{z) (Z = X, Y) is the microscopic density corresponding to the configuration Z: 

pz{z) = J2iS(zi - z) 

• the function w(r) is a short range sigmoid (or step) function such that if r denotes 
the typical radius of the particles, w(r) is close to 1 for r < ar and close to zero 
otherwise. The value of a = 0.3 gives a notion of overlap not too sensitive to small 
atomic displacements. 
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Notice that with this definition q(X, Y) is maximal if X = Y, while it is equal to zero if 
X and Y are uncorrelated. Notice also that one can write q(X,Y) = J2i,j w (\ x i ~ Uj\)i 
a form which is manifestly invariant under permutations of the particles. 

Suppose now to fix a reference configuration Y, chosen with Boltzmann probability 
at temperature T, and consider the conditional probability 

Pq{xlY) = <-"<*W)-«) (6) 

where the constrained partition function is 

Z q {Y) = J dXe- pH ^S(q(X, Y) - q). (7) 

As Y by hypothesis is an equilibrium configuration, it will belong to some quasistate a, 
so, if we choose q as the typical overlap among configurations in this quasistate (with 
probability one almost all configurations have the same overlap), i.e. the Edwards- 
Anderson parameter of the state qEA, we will be able to compute the quasistate averages: 
given two observables A(X) and B(X) we can write: 

Ew a {A) a (B) a = J dY e -^A{Y) J dX e ' mX)6 ^ Y) - q) B(X)(8) 

Notice that, if on the other hand in (|6]) we would choose q as the typical overlap among 
different quasistates, the constraint would be completely irrelevant and we could get the 
Boltzmann average (^). 

Notice that the overlap we consider, is a masked integral of the density-density 
correlation function among the two configurations X and Y. We are interested to study 
the fluctuation of this quantity, which is the following integral of the 4-point function: 

X 4 = P(W - ¥ 2 ) 

= J dx dy dz drw(x — y)w{z — r) 

(((px(x)-p)(py(y)-p))((py^)~p)(Pz(r)-p)} 

-({ P x{x)-p){ P y{y)-p)) x ((px(z)-p)(Pv(r)-p))) (9) 

where we have denoted with the angular brackets the conditional average with the 
distribution and with the bar the average over the canonical distribution of the 
reference configuration Y. 

Within this formalism, the generating functional of the correlation functions is the 
constrained free-energy 

V(q) = -(3/N j dY^- hg(Z g (Y)). (10) 

This function has been computed in various models having a glass transition, 
including mean field disordered models and simple liquids in the HNC approximation, all 
giving consistent results |L7], [18| . The shape of V as a function of q allows to distinguish 
among liquid and glass phase. We show the potential in figure [l] for hard spheres in the 
HNC approximation. 
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Figure 1. The potential for HNC hard spheres at several densities. We show the 
high and low q regions of the potential. The lines joining them are just guides for 
the eyes. Lower curves correspond to higher densities. At low density the potential is 
convex. The appearance of a secondary minimum signals the breaking of ergodicity, 
with exponentially many states. 



In the liquid phase at high temperature, the potential is convex with a unique 
minimum for q = which corresponds to the typical overlap among random liquid 
configurations. Lowering the temperature, the potential looses the convexity, until, 
when T c is reached, it develops a secondary minimum, at a high value of q. The height 
of the secondary minimum with respect to the first one is related to the configurational 
entropy by: V sec — V pr i = TE(T), which vanishes at Tk- A remarkable fact that has 
been often discussed is that while the properties of the low q minimum reflect the 
properties of the full Boltzmann average, the properties of the high q minimum reflect 
the properties of averages in a single ergodic component. 

In the shape of the potential the MC transition appears as a spinodal point and 
as such it has a divergent susceptibility. In fact, general relations in the effective 
potential theory, imply that the susceptibility is given just by the inverse curvature 
of the potential in the minimum, i.e. \i = 1/V"(q) | secondary minimum- This quantity, 
diverges for T — > T c , which, in turn, implies the divergence of the spatial range of the 
correlations. Generically, in all the model studied, the slope of the flex vanishes linearly 
for T — ■> T c , implying Xa{T) — \T — T c |~ 7 with a mean field exponent 7 = 1/2. We 
notice that %4 computed in the primary minimum represents the Boltzmann average of 
the order parameter fluctuations and is completely regular at T c . 

In figure || we see the susceptibility computed with this procedure (circles), which 
shows a divergence at T c . 
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Figure 2. The susceptibility of the metastable states for the p-spin model with p = 3. 
Here T c = 0.612. The low temperature data come from the potential theory. The high 
temperature data from the dynamical equations. 



3. A dynamical approach 

The idea of considering a system coupled with a reference configuration can also be 
used in dynamics to compute the time dependent susceptibility. In this context it is 
convenient to couple with the initial configuration X . Consider a system at equilibrium 
at time zero with respect to the Hamiltonian H(X) which evolves for positive times with 
the modified Hamiltonian 

H tot (X) = H(X)-eq(X,X ). (11) 
For small e, linear response theory at equilibrium implies that 

X4 (t) = (3 ((q(X t , X ) 2 > - (q(X t , X )q(X , X )>) = d M X £ X °» (12) 

The problem of studying the evolution of a system with Hamiltonian (|TT|) can in 
principle be issued within any dynamical approximation scheme (e.g. MCT). 

However for the time being we have only addressed the problem in the context of the 
p-spin model, which, for all the present purposes should capture the essential features 
of the function x^)- Clearly, lacking completely in the model any spatial structure, in 
order to infer from the behavior of Xa something about a correlation length we need to 
resort to eq. ||. 

The p-spin model |T9| describes N interacting variables S\,...,Sn (spins) on 



the sphere J^i.Sf = N, with Hamiltonian H = Y,n<...<i Ji 1 ...i p Si 1 .-.S ip where the 
couplings are random independent Gaussian variables with zero mean and variance J 2 = 
p\ / \ 2N p ~ l ) . The appropriate notion of overlap for this system is q(S, S') = 1/NJ2% SiS'i- 
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For this model it is customary to consider Langevin dynamics, which in our case will 
be performed with Hamiltonian H tot (S) = H(S) — eq(S, S ), where So = S(t = 0) is an 
equilibrium initial condition. 

dSi{t) dH(S(t)) 

— = _ + eS,(0)+*(f) (13) 

where r)i(t) is a white noise with amplitude 2T and fi(t) is a Lagrange multiplier which 
ensures the spherical constraint at all times. 

Using standard manipulations based on the Martin-Siggia-Rose functional integral, 
one can write a self consistent equation for a single spin which, using the notation 
f(q) = l/2q p , reads: 

dS ^ - -fi(t)3(t)+ f ds f"(C(t,s))R(t,s)S(s) 



dt 

+ (3f'(C(t,0))S(0) + eS(0)+Ut) (14) 
where is a colored Gaussian noise with variance 

mC(s)) = f(C(t,8)) + ZT8(t-a) (15) 

where C and R are the correlation and response functions of the system, to be 
determined self-consistently by C(t,s) =< S(t)S(s) >, R(t,s) = The detailed 

derivation of fll4|) is rather standard (see e.g. fl20fl ) and we do not reproduce it here. 

From (|TJ]), taking the correlations with S(s) and £(s) one can derive equations for 
C and R, which read, for t > s 

0C(t - s) -n(t)C(t,s)+ t duf"{C{t,u))R{t,u)C{u,s) 



ds Jo 



+ / f'(C(t,u))R(s,u) + Pf'(C(t,0))C(s,0) + eC(s,0) 

dR{t, s) 







-u(t)R(t,s)+ t duf"(C(t,u))R(t,u)R(u,s). (16) 

OS Js 

Together with the equation specifying the time dependence of /x(£) 

fi(t) = f duf"(C(t,u))R(t,u)C(u,t) 
Jo 

+ f f{C{t, u))R{t, u) + (3f'(C(t, 0))C(t, 0) + eC{t, 0) + T. (17) 
Jo 

Equations (0,0) f° rm a complete set, that can be solved numerically, to derive the 
value of Xi(t) as 

X4 (t) = dC(t,0)/de. (18) 



With a simpl step-by-step integration [21 1 we could reach times of the order of 1000. 



This allowed us, in the case p = 3, to compute the function Xiif) down to temperature 
T = 0.7, compared with a critical temperature Td = 0.612. More sophisticated 
algorithms (see the contribution of Latz to these proceedings) will allow in the next 
future to approach the critical temperature much more. Our results for the function 
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Figure 3. The dynamical susceptibility for the p-spin model (p=3) as a function of 
time for several temperatures the lower is the maximum the higher is the temperature. 
From bottom to top T = 1.0, 0.9, 0.8, 0.7. 

X±(f) are displayed in figure for various temperatures. We see that x nas a maximum 
which becomes higher and higher as the temperature is lowered, and is pushed towards 
larger and larger times. This is the behavior that qualitatively is seen in the numerical 
simulations || [l2| . 

As shown in ref. we define t* as the time at which x* is maximum, we find that 
Xt = Xi(t*) exhibit a divergence at T c , as it is presented in figure ||. Both quantities 
behave as powers of T — T c t* ~ (T — T c )~ a , xt ~ (T — T c )~ 7 . A best fit gives the values 
7 = 0.52 ±0.02, a = 1.1 ±0.1. 

4. Conclusions 

In this paper we have reviewed the analysis of the non linear susceptibility in supercooled 
liquids and glasses that comes from the mean field theory of disordered systems, and 
liquid models in the HNC approximation. 

The theory predicts that while the long time, equilibrium susceptibility remains 
finite and is regular at all temperatures, the finite time susceptibility displays a 
maximum as a function of time which becomes higher and higher and displaced to 
larger and larger times for temperatures close to T c . This behavior is a consequence of 
the critical character of the Mode Coupling like dynamical transition predicted by the 
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ideal theory described in this paper. I real systems one can expect a similar behavion 
but with a round off of the divergence. 
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